//============================================================================
// Name        : TestThreads.cpp
// Author      : Massimo Marchi
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#define _CRT_SECURE_NO_WARNINGS

#include <cmath>
#include <iostream>
#include "Array.h"
#include "fftw++.h"
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>
#include <omp.h>

#include <boost/chrono.hpp>

#include <iomanip>
#include <ctime>
using namespace boost::chrono;

using namespace std;
using namespace Array;
using namespace fftwpp;

using  boost::uniform_real;

template< class Clock >
class timer
{
  typename Clock::time_point start;
public:

  timer() : start( Clock::now() ) {}

  typename Clock::duration elapsed() const
  {
    return Clock::now() - start;
  }

  double seconds() const
  {
    return elapsed().count() * ((double)Clock::period::num/Clock::period::den);
  }
};

boost::mt19937 gen;

const unsigned int NMAX=256;
int main() {
	unsigned int nnx=NMAX;
	unsigned int nny=NMAX;
	unsigned int nnz=NMAX;
	unsigned int nzp=nnz/2+1;

	size_t align=sizeof(Complex);

	array3<double> ro;
	array3<double> ro_out;
	array3<Complex> rok;
	ro.Allocate(nnx,nny,nnz);
	ro_out.Allocate(nnx,nny,nnz);
	rok.Allocate(nnx,nny,nzp,align);
	cout << ro.Nx() << " " << ro.Ny() << " " << ro.Nz() << "\n";
    uniform_real<double> dist(0.0,1.0);

    int op;
#pragma omp parallel
    op=omp_get_num_threads();
/*
    cout << " nthread? ";
    cin >> op;
*/

    fftw_init_threads();
	fftw_plan_with_nthreads(op);
	rcfft3d Forward(nnz,ro,rok);
	crfft3d Backward(nnz,rok,ro);

    for(unsigned int i=0;i<nnx;i++)
    	for(unsigned int j=0;j<nny;j++)
    		for(unsigned int k=0;k<nnz;k++)
    			ro[i][j][k]=static_cast<double> (i);

	cout << "opp " << ro[NMAX/2][NMAX/2][NMAX/3] << endl;

//    cout << " nthreads? " ;
//    cin >> op;

	cout << "lollo " << endl;
	timer<system_clock> t1;
	timer<steady_clock> t2;
	timer<high_resolution_clock> t3;


	cout << "opp " << ro[NMAX/2][NMAX/2][NMAX/3] << endl;

	for(int uo=0;uo<4;uo++){
		cout << op << "\n";


		cout << "nthreads =  " << op << endl;
		cout << "opp " << ro[NMAX/2][NMAX/2][NMAX/3] << endl;

		Forward.fft(ro,rok);

		Backward.fftNormalized(rok,ro_out);

		cout << "ipp " << ro_out[NMAX/2][NMAX/2][NMAX/3]  << endl;

	}
	std::cout << "system_clock-----------: "
			<< t1.seconds() << " seconds\n";
	std::cout << "steady_clock--------: "
			<< t2.seconds() << " seconds\n";
	std::cout << "high_resolution_clock--: "
			<< t3.seconds() << " seconds\n";

/*
#pragma omp parallel num_threads(3)
    {
	#pragma omp parallel for
    	for(int n=0; n<size; ++n){
    		sinTable[n] = std::sin(2 * M_PI * n / size);
    	}

    	// the table is now initialized
    	cout << "!!!Hello World!!!" << endl; // prints !!!Hello World!!!
#pragma omp master
    	{
    		int nthreads=9;
    		cout << "There are " << nthreads << " threads" << '\n';
    	}
    }
*/
	return 0;
}
